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We consider polarons in models of coupled electronic and vibrational degrees of freedom, in the 
presence of a soft nonlinear inter-particle potential (Morse potential). In particular, we focus on a 
a bound state of a polaron with a breather, a so-called "polarobreather" . We analyze the existence 
of this branch based on frequency resonance conditions and illustrate its stability using Floquet 
spectrum techniques. Multi-site solutions of this type are also obtained both in the stationary case 
(two-site polarons) and in the breathing case (two-site polarobreathers). We also obtain a different 
branch of solutions, namely a polaronic nanopteron. 
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I. INTRODUCTION 

In the past few years, the study of discrete nonlinear 
systems has received a lot of attention, chiefly due to the 
increasing number of pertinent physical applications 0. 
Among the various areas of recent interest, one can high- 
light coupled waveguide arrays in nonlinear optics 
Bose-Einstein condensates (BECs) trapped in deep opti- 
cal lattices (OLs) [jj in atomic physics, coupled cantilever 
systems in nano-mechanics |^| , the local denaturation of 
the DNA double strand in biophysics [f| , stellar dynamics 
in astrophysics 6], and so on. 

In the study of such discrete nonlinear systems, an im- 
portant concept is the possibility of intrinsic localization 
due to interacting degrees of freedom, even in a linear 
regime, as proposed by Holstein This type of sys- 

tems has been shown to sustain single- and multi-humped 
polaronic and excitonic solutions |£g. The same notion 
of effective nonlinearity in the setting of coupled exci- 
tonic and vibrational degrees of freedom was central to 
Davydov's suggestion of solitonic excitations arising in 
biomolecules [lol Hi) . More recently, a new aspect has 
been added to this type of problems 0, 0] , by consid- 
ering the interplay of the linear self-trapping with a soft 
nonlinear potential, such as a Morse potential. This has 
been proposed as a more general model relevant to many 
soft matter applications. 

Our purpose in the present paper is to examine this 
type of model and to provide a number of insights re- 
garding its solutions and their stability. We start by es- 
tablishing that the main stationary solutions obtained 
in the earlier work of [T^ (namely, the single-site po- 
larons) are dynamically stable within their region of ex- 
istence. We then turn our attention to a different class 
of solutions which are genuinely "breathing" in time and 
which will, thus, be termed polarobreathers. Such so- 
lutions were first obtained and discussed in the context 
of the Holstein model in [Til ITU , from where this term 



is coined. Here, we will obtain (and justify, based on 
resonance conditions) the domain of existence of such so- 
lutions. We will also systematically investigate their sta- 
bility, by performing the corresponding Floquet spectral 
analysis, and illustrate the existence of intervals of linear 
stability of these solutions. Our study of this new class of 
solutions will reveal yet another family of solutions con- 
sisting of a static polaron combined with a linear mode 
to give a form of polaronic nanopteron. Since the branch 
of stationary polarons terminates at a critical value of 
the coupling as illustrated in the earlier work of fl2^ , we 
will examine the dynamics beyond this termination point 
with a single hump initial condition and will show that 
the solution develops large breathing fluctuations in its 
local (i.e., central site) energy in time due to an inter- 
play between the central site of the localized profile and 
its neighbors. We will also highlight multi-site branches 
such as the two-site stationary solutions (polarons) and 
the two-site polarobreathers and numerically reveal their 
instability. 

Our presentation will be structured as follows. In sec- 
tion II, we will present the model and numerical meth- 
ods (the details of which are relegated to a technical ap- 
pendix). Subsequently, in section III, we will present and 
discuss our numerical results. Finally, in section IV, we 
will summarize our conclusions and pose questions of in- 
terest for future studies. 



II. THE MODEL 

We consider the coupled charge/excitation-lattice 
model introduced in describing the competition be- 
tween linear polaronic self trapping and self-focusing ef- 
fects of a soft nonlinear potential; in dimcnsionlcss form, 
it can be expressed as follows: 

ii> n = -J(*„+i + *„-i) - FA, (1) 
u n = -V'(u n ) + xl*n| 2 + k(u n+1 - 2u n + Un-x), (2) 
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where dots denote time derivatives, and the lattice index 
n runs from 1 to N (the total number of lattice sites). In 
Eqs. GJ-0, ^n(t) represents the "electronic" degrees of 
freedom, u n corresponds to the lattice displacements (i.e., 
"vibrational" degrees of freedom), while the parameters 
J, k and x denote, respectively, the transfer integral, the 
lattice spring constant and the coupling constant between 
the interacting fields. Finally, V(u n ) is an anharmonic 
on-site potential, which, similarly to 0], is assumed to 
have the form of a Morse potential: 



V(u) = -{eM-u)-l} 2 . 



(3) 



The dynamical equations iJTJ-© arise from the Hamil- 
tonian: 



1 k 

-ii 2 n + V{u n ) + -(u n - u n+1 ) 2 



H = J2 

n 

- E [X(\*n\ 2 u n ) + J(*nK + l + C -C-)] (4) 



We will seek standing-wave, as well as genuinely time- 
periodic solutions of Eqs. CJ-©- To factor out the phase 
invariance of Eq. (JIJ, and clearly separate between the 
two classes of solutions, we introduce the transformation: 



H„(t) = * n (t)e-*-*. (5) 

Then, the dynamical equations take the following form: 

lE n - cj c S„ + J(H„ + i + H n _i) + X u n^n = 0, (6) 
u n + V'(u n ) - x|S„| 2 - k(u n+ i - 2u n + u n -i) = 0. (7) 

Stationary solutions now correspond to ones where all 
time derivatives in Eqs. ©-(0 are se t to zero. On the 
other hand, time periodic solutions are characterized by 
a frequency of oscillation lu^ (or period Tb = 2ir/u>t,) 
in the time dependence of both the electronic wavefunc- 
tion S„ and of the lattice displacements u n . Our aim 
is to study the existence and stability of such stationary 
and breathing polaron solutions; the latter, adopting the 
terminology of (l4Lll5j . will be henceforth called "polaro- 
breathers" . In the results to be presented below, we have 
fixed the values of the transfer integral and lattice spring 
constant, namely J = 0.005 and k — 0.13, and let the 
coupling constant x an d the polarobreather frequency 
Wb vary. Nevertheless, in our study, we have also con- 
sidered other parameter values (such as, e.g., k = 0.065, 
k = 0.26, J = 0.0025 and J = 0.01) and we have obtained 
qualitatively similar results to those reported below. 

In order to find either stationary or time-periodic solu- 
tions of the dynamical equations ©-Q), we use methods 
based on the anti-continuous limit [T(j,ll2l- In particular, 
we first obtain a nontrivial steady state uq -for station- 
ary solutions- or an orbit of frequency LOh for an isolated 
oscillator uo(t)-for time periodic solutions. Then, the 
solution at the anticontinuous limit (k = J = x = 0) 



is S„ = u n = 0, except at n = 0, which is set to u$. 
The coupling constants are subsequently varied through 
a path-following (Newton-Raphson) method. In the fol- 
lowing, we first vary x up to 0.02, then J up to 0.005 and, 
finally, k up to 0.13. Once J and k are fixed, x can be 
freely varied. The implementation of the anticontinuous 
limit has been performed in real space (using a shoot- 
ing method) and in Fourier space. Note that the first 
method is less accurate than the second one, whereas the 
latter is rather time consuming, particularly in the case 
of long chains. We have made use of both methods in 
the calculations presented herein. 

Once the numerically exact [up to the prescribed ac- 
curacy O(10 -7 )] solutions are obtained, linear stability 
analysis is performed to examine the dynamical stabil- 
ity of the solutions. In the case of stationary solutions, 
this is implemented by imposing normal mode pertur- 
bations (see, e.g., earlier work in 0, 0, 0]) and ob- 
taining the corresponding eigenvalues of the ensuing lin- 
ear matrix problem. On the other hand, in the case 
of (time-periodic) polarobreather solutions, we impose 
time-dependent perturbations to both the electronic and 
vibrational fields and solve the resulting (linear) differen- 
tial equations for the perturbations from t — to t = Tt,. 
Then, obtaining the Floquet (or monodromy) matrix, 
which relates the perturbation vector at t = to that 
at t = Tb, we compute the Floquet multipliers of the pe- 
riodic solution. The Floquet multipliers imply stability 
(instability), if they do (do not) appear only on the unit 
circle. The details of the numerical methods used in this 
work are discussed in the Appendix. 



III. NUMERICAL RESULTS 

In this section we illustrate the numerical results ob- 
tained for the polarobreathers. Recall that in all calcu- 
lations we have fixed J = 0.005 and k = 0.13, x being a 
free parameter. 

We first consider the stability of static polarons. Such 
polarons, localized on a single lattice site, are stable. 
Considering perturbations of these structures, namely 

n,0 



«n(0) 



l(0)+e<5„,o andi?e(*„(0)) = #e(*° (0))+e5, 



(where the zero superscript denotes an exact solution), it 
is found that the polaron evolves into a vibrating state, 
which is a normal mode of the system. This result is 
clearly illustrated in Fig. ^(here we have used e = 0.001 
and x = 0-5); note that in the relevant Fourier spectrum 
(right panels of the same figure for the vibrational (top) 
and electronic (bottom) field at the central site), there 
appears a single excited frequency. 

Apart from static polarons, there exist stable, time- 
periodic and spatially localized solutions within the con- 
text of the model, namely polarobreathers. Fig. [21 shows 
an example of these solutions along with the relevant Flo- 
quet spectrum. Additionally, Fig. [21 shows the Fourier 
spectrum of this solution, whose nonlinear character im- 
plies the existence of integer multiples of its frequency. 
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FIG. 1: Left panel: Evolution of the displacement of the cen- 
tral particle of a slightly perturbed static polaron with \ = 0.5 
Right panels: The Fourier spectrum of the corresponding time 
series is shown (the top part corresponds to the lattice vari- 
able and the bottom part to the electronic one). The existence 
of only one frequency indicates the simple periodic nature of 
the solution. 
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FIG. 2: Left panels: Profiles of the lattice (u n , top) and elec- 
tronic component ( | 2 , bottom) of a polarobreather with 
ujb = 0.8 and x — 0.15. Right panel: Floquet spectrum of 
this solution. 
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FIG. 3: Fourier spectrum of the polarobreather of Fig. |5] The 
top and bottom panels correspond to the lattice and electronic 
variables respectively. 
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FIG. 4: Time evolution of the lattice (left panel) and elec- 
tronic coordinates (right panel) of a polarobreather with 
u>b = 0.8 and \ — 0.15. Notice that the oscillation period 
of the two functions is different. 



The time evolution of the lattice (uo(t)) and electronic 
(^o(i)) variables of the central particle is displayed in 
Fig. 01 Notice the different period in the oscillation of 
the u and the \& fields, due to the additional inclusion, in 
the latter, of the oscillation with frequency u> e . Polaro- 
breathers exist as long as the conditions of the MacKay- 
Aubry's theorem [16j are fulfilled, i.e., none of the har- 
monics of the polarobreather frequency resonate with the 
linear modes. Based on this criterion, an analysis of the 
linear modes can provide the range of existence of po- 
larobreathers of a given frequency. Fig. shows the 
real part of the normal mode frequencies and indicates 
the existence of a continuum band of (extended) linear 
modes besides several localized modes. Since we are deal- 
ing with polarobreathers in a soft potential, they must 
stem from the linear mode at the bottom of the band of 
extended modes. This mode becomes localized (i.e., bi- 
furcates from the continuous spectrum band) at the value 
of x ~ 0.033 and its location varies with x- As a result, 
the polarobreather frequency must be smaller than this 
localized-mode frequency. This fact is indicated in the 
right panel of Fig. [S] It can be the observed that the bi- 
furcation points correspond to the values of x f° r which 
the frequency of the localized mode coincides with the 
frequency of the polarobreather branch. 

The existence of polarobreathers is also limited by the 



second-harmonic resonances. In the case of finite lat- 
tices, "phantom" solutions can appear due to the exis- 
tence of gaps in the extended modes band These 
solutions consist of a core that vibrates with the funda- 
mental frequency uj^ (or may also be static), whereas the 
tails correspond to a phonon vibrating with frequency 
2wb- These solutions are an artifact of the finiteness of 
the lattice and, moreover, they are usually unstable. An 
example of a phantom polarobreather is shown in Fig. 
(for = 0.6 and x = 0.15). 
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FIG. 5: The left panel shows the real part of the linear mode 
spectrum. The right panel shows the dependence of the po- 
larobreather energy as a function of x f° r different values of 
the frequency (from left to right): cu h = 0.9, 0.84, 0.8, 0.7. The 
line at the bottom corresponds to the static polaron. 
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FIG. 6: Left panels: The profiles of the lattice (top) and the 
electronic (bottom) components of a phantom polarobreather 
are shown for ujb = 0.6 and \ = 0.15. Right panel: Floquet 
spectrum of this solution. 
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FIG. 9: Floquet spectra for a two-site polarobreather, subject 
to a harmonic (left panel, \ — 0-5 ) or an oscillatory (right 
panel, x — 0.15) instability. 
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FIG. 7: The argument (left) and modulus (right) of the Flo- 
quet eigenvalues as a function of the coupling constant \ f° r 
one-site polarons. 




An analysis of the stability of polarobreathers is shown 
in Fig. [7| It is observed that, for the considered branch, 
the only kind of instabilities that arise are oscillatory ones 
(stemming from complex Floquet multipliers). The small 
growth rate of the ensuing instabilities, however, ren- 
ders the polarobreathers relatively long-lived metastable 
structures when they are not linearly stable. 

Apart from one-site (site-centered) polarons, we have 
also considered two-site (bond-centered) static polarons 
and polarobreathers. Whereas one-site polarons are sta- 
ble, two-site ones are unstable. An example is shown in 
Fig. [8] where the real and imaginary parts of the linear 
mode spectrum for two-site polarons are illustrated. 

Similarly to their static counterparts, two-site polaro- 
breathers are unstable. This instability can be either 
exponential (i.e., the Floquet eigenvalues responsible for 




FIG. 8: The real (left) and imaginary (right) part of the linear 
mode spectrum as a function of the coupling constant x f° r 
two-site polarons. 



FIG. 10: Domain of existence of one-site (black region) and 
two-site (black and gray region) polarobreathers in the x~ w b 
parameter plane. 

the instability have zero phase) for small Xi or oscillatory 
(i.e., the Floquet eigenvalues responsible for the instabil- 
ity have phases different from zero or n). In Fig. [5J 
examples of the Floquet spectra pertaining to a regular 
exponential (for x — 0.5) and an oscillatory instability 
(for x — 0.15) are shown. 

Two-site polarobreathers can be continued as a func- 
tion of x f° r J and k fixed. As in the one-site case, the 
branches for a fixed frequency merge with a static solu- 
tion branch. The domains of existence of the one-site and 
two-site polarobreathers for k = 0.13 and J = 0.005 are 
shown in Fig. ^| It is interesting to note that, similarly 
to what happens for their static counterparts [l^, one- 
site polarobreathers have a narrower domain of existence 
(gray region in Fig. HOfl than their two-site counterparts 
(gray and black region in Fig. I1U|) . Given that the termi- 
nation of such branches occurs upon their collision with 
the stationary branch of solutions, this trait is rather nat- 
ural to expect in the present setting. A wider domain of 
existence for two-site breathers was also observed in a 
model with competing attractive and repulsive interac- 
tion m. 

In order to examine the dynamical evolution of the 
instability of two-site polarons, we again consider the 
perturbation u n (0) = it°(0) + £# n ,o and Re(9 n (Q)) = 
Re(^(0)) + eS n ,o, with the value e = 0.001 pertaining 
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FIG. 11: Evolution of the central particles displacements and 
energy densities, respectively, for a perturbed two-site polaron 
with x — 0.35 (left panel) and a two-site polarobreather with 
X = 0.2 (right panel). 




FIG. 12: (Color online) Left panel: Time evolution of the 
displacement of the central particle at \ = 0.72 using as initial 
condition a static polaron at \ — 0.718. Right panel: the 
energy density plot. 



to a static polaron. It can be observed that, as illustrated 
in Fig. 1111 after a transient, the polaron evolves into an 
anti-phase vibrating state; note that if the two central 
particles were perturbed, the observed behaviour would 
be similar, but the transient would be longer. If the same 
perturbation is applied to a two-site polarobreather, the 
behavior is somewhat reminiscent of that observed for a 
static polaron. This is illustrated in Fig. ^2 by display- 
ing the temporal evolution of the energy density at the 
two principal sites of the solution. The energy density is 
defined as: 



1 k 

-U% + V(u n ) + - [(ti„ - U n+ if + (u n - U„-i) 2 ] 

X |*„| V - ^ n (K+i + K-i) + cc.) (8) 



IV. CONCLUSIONS 

In this work, we have discussed the existence and 
stability of localized solutions in a generalized Holstcin 
model with a soft (Morse-type) inter-particle potential. 
We have shown that in addition to the standard station- 
ary polaronic solutions within this model discussed else- 
where |l2j , there exists also very robust breathing behav- 
ior in the form of a time-periodic type of state: the po- 
larobreather solutions. Using Floquet analysis, we have 
illustrated the existence of parametric intervals of stabil- 
ity of the latter solutions, and we have also elucidated 
their domain of existence based on resonance conditions. 
In this respect, we have shown that the polarobreathers 
cease to exist as a result of their "collision" with the stan- 
dard polaronic branch, occurring at frequency-dependent 
values of the relevant coupling parameter. We have also 
found that "phantom" versions of such polarobreathers 
exist in finite lattices but are usually unstable. We can 
view these dressed polaron solutions as "hot" polarons 
(or excitons in related exciton-vibron systems). However, 
the polarobreather is dressed with a coherent phonon 
wavepacket (a breather), which binds with the polaron 
as a composite local excitation. We have also obtained 
multi-site polarons (such as the two-site polarons ex- 
amined herein) and multi-site polarobreathers (such as 
the two-site polarobreathers, also studied herein), and 
numerically demonstrated their instability. Finally, we 
have shown that the unstable polaronic dynamics (for 
parameter values beyond the range of existence of their 
corresponding stationary branch) may lead to energy re- 
distribution and eventual oscillatory dynamics for longer 
timescales. 

It is interesting to examine higher-dimensional gener- 
alizations of the model and the possibility for obtaining 
localized solitary wave or vortex-like solutions in such 
higher-dimensional discrete contexts. Such studies are 
currently in progress and will be reported elsewhere. It 
would also be important to assess the significance of 
quantum effects in both the electronic and lattice de- 
grees of freedom. Initial results [2l| suggest that bound 
states of breathers and polarons (or excitons) may indeed 
persist in appropriate regimes. 



Another interesting fact is that the static polaron 
branch does not exist for \ > Xc ~ 0.718. In particu- 
lar, if the solution for a polaron close to \ c is used as 
initial condition for \ > Xc), e.g., for \ — 0.72, we have 
found the following type of evolution (see figure lT2"|) : first, 
for some finite time, a continuous increase of the "local" 
energy at the central site of the solution occurs (coupled 
to a corresponding decrease in the energy of the imme- 
diate neighboring sites). Then, the system "relaxes" so 
that a quasi-breathing evolution of the energy density is 
observed; however, its strong localization at the central 
site of the lattice is essentially preserved for long times. 



V. APPENDIX: DETAILS OF NUMERICAL 
METHODS 

A. Real space methods 

In order to analyze Eqs. ©-(0 using the real space 
method, it is convenient to separate the electronic wave- 
function into its real and imaginary part, 5 n (t) — 
4>n(t) + i(p n (t), and also define the oscillation velocity 
v n (t) — u n {t). In this way, we end up with a system of 
AN equations (N is the number of particles of the sys- 
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tern): 

4>n = -(-U e +XUn)<Pn — J(<Pn+l+<Pn-l), 
fin = (—W e + X u n)<t>n ~ J {4>n+l + 4>n-l), 

V n = -V'(u n )+ X (<t> 2 n + V 2 n) 

+ k(u n+ i - 2u n + u n -x). (9) 

Then, if we define X(t) = {u n (t),v n (t),<p n (t),ip n }, there 
exists a map T which relates X(0) and X(Tb) as X(J\,) = 
TX(0). The shooting (real space) method consists in 
finding zeros of the map F = TX(0) - X(T b ) [13 (f or 
a detailed explanation of the shooting method, see |22j). 
The tangent map needed for the application of the New- 
ton method, J = dT, must be calculated numerically. 

However, in this scheme the electronic norm is not con- 
served and the electronic frequency u e must be known. 
For this reason, an equation is added to |(3J|, 

+ -i = o, (io) 

n 

so that the norm is conserved, and a new variable lo c is 
added to the vector X(t). 

B. Fourier space methods 

The Fourier space methods (see details in [23, ) are 
based on the fact that both the electronic wavefunction 
and the lattice displacements are periodic with period 
Tb. Thus, they can be expressed in terms of a truncated 
Fourier series expansion: 

k m 

k=-k m 
k m 

«n(*) = b n C M ik ^bt). (11) 

k=-k m 

Then, the dynamical equations ©-0 are reduced to 
a set of (2N+1) x (2k m + 1) algebraic equations (we have 
also included the norm conservation) , where the variables 
are Z = {a k ,b k ,u e }-- 

F% n = -uj h ka k n - Lo c a k n + T k (V'(u n )) + x^ k (u n ^n) 

+ J(a n - 1 +b k n+1 ) = 0, 

F k ;7l ee -u;iebt+T k {V\u n ))-xF k {\1 n \ 2 ) 

- k(b k _ 1 -2b k +b k n+1 ) = 

F N , n EE C£J2J2 a n a n)- 1 = - ( 12 ) 

n k k' 

In the above equations, T k denotes the Discrete Fourier 
Transform: 

k m 

F k (u) = ^2 exp(ifcwbtj), (13) 

j=-k m 



where tj is a sample of times that must be equally spaced 
according to 

2nj 

tj = 77yi \ 7\ i 3 = ~k m , . . . , +fc m , (14) 

U\,\£K m ~t~ i-j 

and u(tj) is calculated from the Fourier coefficients b k by 
means of the Inverse Discrete Fourier Transform: 

k m 

u(tj) = b k exp(ikoj b tj). (15) 

k=-k m 

The Newton operator or Jacobian J ee OF (with F = 
{Fp ; Fg ; Fjy} can be calculated analytically. One of 
the main disadvantages of this method is that the Ja- 
cobian is singular. Polarobreathers must be calculated 
using singular value decomposition |25| . However, this 
method allows the calculation of non-time-reversible so- 
lutions. 



C. Normal modes and linear stability analysis 

Here, we present the normal modes and stability equa- 
tions, which were used in the results presented in the 
main text. 

1. Stationary Solution Stability: Normal modes 

In order to calculate the normal modes we first need 
to calculate a stationary polaron from Eqs. Q-© from 
the following conditions |12| : 

*n(*) = ip n exp(-iuj e t), 

u n {t)=u n (t) = Q. (16) 

Then, the polaron is determined by the electronic wave- 
function xp n and the lattice displacements y n . 

The normal modes are introduced as perturbations to 
a stationary polaron [H, 0, 0] : 

%n(t) = bPn + exp(-ZCJ c i), 
Un(t) = y n + fri(*)> 

U n (t) = y n + 7T„(t) = 7T n (t), (17) 

and the dynamical equations for the normal modes are: 

£,n ^n 7 

tin = -V"{y n )£ n + xVVi^n + <), 

+ k(£ n+1 - 2£„ + £ n _i). (18) 

Normal mode frequencies can be calculated using the fol- 
lowing relations: 

£n(t) = o, n exp(iwt) + b n exp(— icu*t) 
£n(t) = c„ exp(jwi) + c* exp(-ioj*t) 
ir n (t) = d n exp(iwi) + d* exp(— iu*t). (19) 
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Thus, the normal modes are determined by the following 
linear system: 

uja n = u e a„ + x{%lnan + IpnCn) + J{a n +i + «n-l), 

LUC n = -id n , 

uid n = i[V"(y n )c n - xi>n(a n + b* n ) 

+ kiCn+i - 2Cn + Cn_i)]. (20) 

As the normal mode operator is non-Hermitian, the 
frequencies are in general complex numbers. If non-real 
frequencies exist in the spectrum, then the stationary 
polaron is unstable. 

2. Polarobreather Stability: Floquet Spectrum 

Linear stability can be studied via a Floquet analysis. 
This analysis is performed by linearizing the dynamical 
equations around a polarobreather. Then, the analysis 
performed for the normal modes cannot be extended to 
the stability, as the electronic wavefunction does not pos- 
sess gauge invariance. Thus, we introduce a perturbation 



n = {a n , f3 n ,£,n,Tr n } to Eqs. 13, which has the form: 

</>n = 4>n + OZ n , <p n = <Pn + fin, U n = U° + £„, 7T„ = £„, 

(21) 

where the superscript zero denotes the polarobreather so- 
lution. Then, to leading-order approximation (i.e., keep- 
ing only the linear terms), the stability equations are: 

&n = -(-W e + XU n )Pn ~ X¥>n£« _ J (<*n+l + OL n -l), 
$n = (-W e + XU n )a n + X<t>n£,n + J (fin+1 + Pn-l), 

ir n = 2 X (^a n + <p° n /3 n ) - V"(u° n ) 

- k(£ n+1 -2£ n +^ 1 ). (22) 

The Floquet operator (or monodromy matrix) M. relates 
the vector 51 at t — and t — Tb as follows: 

Q(T b ) = MQ(0). (23) 

A polarobreather is stable if all the eigenvalues of the 
monodromy matrix lie on the unit circle [lflj . 
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